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We consider possible exotic ground states of quantum spin ice as realized in rare earth pyrochlores. 
Prior work in Ref ll5l introduced a gauge mean field theory (gMFT) to treat spin or pseudospin Hamil- 
tonians for such systems, reformulated as a problem of bosonic spinons coupled to a U(l) gauge field. 
We extend gMFT to treat the most general, nearest neighbor exchange Hamiltonian, which contains 
a further exchange interaction, not considered previously. This term leads to interactions between 
spinons, and requires a significant extension of gMFT, which we provide. As an application, we 
focus especially on the non-Kramers materials Pr2TM207 (TM =Sn, Zr, Hf, and Ir), for which the 
additional term is especially important, but for which an Ising-planar exchange coupling discussed 
previously is forbidden by time-reversal symmetry. In this case, when the planar XY exchange is 
unfrustrated, we perform a full analysis and find three quantum ground states: a U(l) quantum spin 
liquid (QSL) , an antiferro-quadrupolar ordered state and a non-coplanar ferro-quadrupolar ordered 
one. We also consider the case of frustrated XY exchange, and find that it favors a rr-flux QSL, 
with an emergent line degeneracy of low energy spinon excitations. This feature greatly enhances 
the stability of the QSL with respect to classical ordering. 



I. INTRODUCTION 

The quest for quantum spin liquid (QSL) ground 
states, exotic phases of matter with emergent gauge 
structure and quasiparticles carrying fractional quan- 
tum numbers- 1 , is an ongoing endeavour in condensed- 
matter physics. Well studied candidates include some 
two-dimensional organic crystals^ and some inorganic 
kagome systems such as herbertsmithite^. Among three- 
dimensional materials, experimental candidates include 
several magnetic pyrochlore oxides^ and hyperkagome- 
lattice magnets^. Classical spin liquids have been re- 
alized in the spin ices^, in which the spins reside on 
a pyrochlore lattice and interact via a dominant clas- 
sical Ising coupling. It has been shown theoretically that 
a weak quantum-mechanical perturbation does not pro- 
duce long-range order in the ground stated Instead, it 
lifts the macroscopic degeneracy of the spin-ice manifold, 
leaving gapless "photon" excitations, describable by an 
emergent U(l) gauge field. The photon exists in a so- 
called Coulomb phase or U(l) spin liquid phase, which is 
stable to all weak perturbations, at zero temperature. 

To describe the low-energy physics of magnetic py- 
rochlore oxides associated with local magnetic doublets of 
rare-earth ions, a minimal pseudospin- 1/2 model can be 
introduced on symmetry grounds^ (see Eq. ([I])). It has 
also been derived micropco pically using superexchange 
theory for various materials^tL^. This model success- 
fully explains spin correlations experimentally observed 
in Yb2Ti207 (Refs. 8 9). As can be seen from the gen- 
eral form of Eq. (IJ, these comparisons between theory 
and experiment also reveal that putative continuous ro- 
tational symmetry of the pseudospins is broken by a sig- 
nificant level of magnetic anisotropy. Moreover, at least 
for Yb 2 Ti 2 07 and possibly for other materials, the Ising 
interaction remains dominant, in which case the physics 



is that of a quantum variant of spin ice^. At a phe- 
nomenological level, recent experimental findings suggest 
the relevance of the Coulomb phase physics in real rare- 
earth magnetic pyrochlore oxides™^. 

Based on this observation, detailed analyses of the non- 
perturbative stability of the Coulomb phase and the pos- 
sible existence of other phases and phase transitions are 
called for. It must be noted that this is a very com- 
plex problem; the general pseudospin Hamiltonian in 
Eq. contains four exchange constants: the Ising ex- 
change J zz , and three "quantum" terms J±, J z ±, and 
J±±. Assuming we start from the classically frustrated 
spin ice case J zz > 0, one then has three dimension- 
less couplings J±/J zz , J z ±/J zz and J±±/J zz , forming a 
three-dimensional phase space even at zero temperature. 
The development of a comprehensive theory of this full 
3d phase space is a challenging task. 

A method for analysis of this problem was developed 
in Ref. 1 15, based on a gauge theory reformulation of the 
problem on a dual diamond lattice. There, the original 
Hamiltonian was re-expressed as a problem of bosonic 
spinons hopping in the background of a fluctuating com- 
pact U(l) gauge field. This problem was in turn subse- 
quently approximated using a mean-field theory. In that 
work, this gauge Mean Field Theory (gMFT) was applied 
to the corner of the phase diagram approximately appro- 
priate to Yb2Ti207, with, in our (and their) notation, 
J±± = 0, and J± > 0. Both the expected U(l) QSL 
phase and an additional exotic state, a Coulomb ferro- 
magnet, were found, though somewhat limited in their 
domain of stability. 

Here we extend the theoretical formalism to allow to 
fully treat the most generic nearest-neighbor pseudospin- 
1/2 Hamiltonian, i.e. the fully general form of Eq. ([I]). 
This requires some significant technical extensions to the 
analysis in Ref. 1151 In particular, the term J±± induces 
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interactions amongst the spinons, which may induce pair- 
ing and other effects. Furthermore, in the case J± < 0, 
a non-zero average gauge flux is present and complicates 
the dispersion of the spinons. Surmounting these tech- 
nical obstacles, we then apply the extended method to 
the particular case where the lowest crystal-field levels 
of the rare-earth ions are given by non-Kramers mag- 
netic doublets with integer spins. This has a direct 
relevance to PT 2T M2O 7 with tran sition-m etal elements 
TM = Sri™, Zr™l Hf, and I j 16 ' 19 ' 20 '. One key re- 
sult of this analysis is that the U(l) QSL phase is much 
more stable than in the case of a Kramers (half-integer 
spin) doublet, because of the absence of the coupling be- 
tween the Ising and planar components of the rare-earth 
moments (Sec. III). Moreover, the QSL becomes particu- 



larly robust when the U(l)-symmetric planar pseudospin- 
exchange interaction is frustrated (J± < 0). 

The remainder of the paper is organized as follows. In 
Sec. [TlJ we introduce the most generic nearest-neighbor 
pseudospin-1/2 Hamiltonian in rare earth pyrochlores 
and reformulate it as a problem of bosonic spinons cou- 
pled to a U(l) gauge field. Focusing on non-Kramers 
doublet, we show in Sec. Ill a mean-field analysis of the 
gauge theory, following and extending Ref. UM Within 
mean-field analysis, we also discuss that the different 
flux patterns of gauge fields are favored in the presence 
of (un)frustrated planar XY exchange cases, based on 
a perturbation argument. In Sec. |IV[ unfrustrated pla- 
nar XY exchange case (FM case) is studied as follow- 
ing orders: In Sec. IV A| and Sec. |IVB we discuss the 
possible order parameters for spinons and the details of 
self-consistent equations. Then gMFT phase diagram is 
shown in Sec |IV C| with the characteristics of their phase 
transitions described in Sec. |IVD| In Sec [Vj we discuss 
the enhancement of the stability of the QSL for the frus- 
trated planar XY-exchange case (AF case). We finally 
conclude in Sec. | VI| with a summary of results, a discus- 
sion of relevant experiments and open questions. 



II. COMPACT ABELIAN HIGGS MODEL FOR 
PSEUDOSPIN-1/2 QUANTUM SPIN ICE MODEL 

A. Nearest-neighbor pseudospin-1/2 quantum spin 
ice model 
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TABLE I: Local coordinate frames for the four sublattices on 
the pyrochlore lattice. 



a (pseudo-)spin 1/2 operator acting within the Hilbert 
space of the local doublet of the rare-earth ion located at 
that site. Its physical meaning is discussed further below. 
The matrix 7^ = l 7 e ±l27T / 3 (and £ = —7*) depends on 
the bond direction between the neighboring sites i and j 
as seen in Fig JT] 
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Here we defined the four vectors, e p (^ = 0,1,2,3), con- 
necting a site on the diamond sublattice A to its neigh- 
bors. In our coordinates, where the conventional cubic 
unit cell is taken of length 1, these are related to the local 
z qua ntiza tion (Ising) axis for the pseudospins given in 
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Now we review the physical meaning of the pseu- 
dospins. The z component Sf is directly proportional to 
the physical magnetic moment along the local Ising (111) 
axis. The interpretation of the "in-plane" pseudospin 
operators, Sf,Sf, differs, however, for the Kramers and 
non-Kramers cases. In the former, this is indeed pro- 
portional to the magnetic dipole moment normal to the 
local (111) axis. Because in the non-Kramers case, the 
in-plane components of the pseudospin are time-reversal 
invariant, they must be identified not with the magnetic 
dipole moment but for instan ce the quadrupole moment 
Sf ~ { Jf , Jf } in the Pr cases!™. 

This is summarized formally as follows. Defining lo- 

and Zi 



cal coordinate axes Xi, y i 
magnetic dipole moment is given as 



as in Table II A the 



We start with the generic nea rest- neighbor pseudospin- 
1/2 model for quantum spin ic d 8 * 11 * ^, 

H = Y,l J *» S " S 'j ~ J ±( s * +S 7 + S * rS i~) 

m 

+J ±± (7i i S+Si-+ 7 £S i -ST) 

+J Z± (SKC^S+ + QSj) + (C 4J S+ + £.Sr)S<)]. 

(1) 

Here, Y^Uj) indicates the summation over all the nearest- 
neighbor sites i and j on the pyrochlore lattice. Si is 



(Mi) = //.„/'/;; S; x, + (SDi/i) + g z u. B {Sf)Zi, (3) 

where for the non-Kramers magnetic doublets which are 
of our particular interest in this paper, g xy = 0. For this 
case (and specifically for Pr 3+ ) the quadupole moment 
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FIG. 1: Mapping onto U(l) gauge theory : spinons are denned 
on diamond lattice sites (black sphere) and gauge fields live 
on their links (black solid line). J± (J±±) shows quadratic 
(quartic) spinon hopping process in Eq.(|ll[). 



B. Interacting bosonic spinons coupled to compact 
U(l) gauge fields on a diamond lattice 

In this section, we exactly re-express the Hamiltonian 
Eq. as a U(l) gauge theory on the dual diamond 
lattice. The formulation is identical, apart from some 
minor notational changes, to the one in Ref. 1151 except 
that we now include the J±± term. 

First, the Ising interaction term is simpl y exp ressed as 
Jzz J2 r Qr w i tn a g au g e ( "electric" ) charge^Zl 

3 

Qr = Vr ^ (5) 

where the coordinate r labels the centers of the py- 



rochlore tetrahedra, which form the "dual" diamond lat- 
tice. We furthermore defined rj r = 1 and -1, when r is 
on the A and B diamond sublattice, respectively, (see 
FigJT]) Note that the electric charge Q r takes integer val- 
ues. Then, the corresponding "electric field" can be taken 
as a directed link variable, 

i'r,r+j )r e A , = Vr^l+ Vrefi /2 = 'fc S r,r+»7 I .e (J ■ (6) 

With these definitions, Eq. ^ may be regarded as the 
lattice analog of Gauss' law. Maintaining the constraint 
m Eq. (§, one then represents the planar components 
of the pseudospin operator as 

S r + +e (1 /2 = ^r+e^r+e, for r G A, (7) 

where $ r and &l are annihilation and creation operators 
of bosonic spinons that decrease and increase the "electric 
charge" Q r , respectively, 

[$r,Or]-<f>r, H, Qr] = (8) 

Note that conceptually plays the role of the expo- 
nential of the gauge vector potential A, i.e. an element 
of the gauge group, which creates or annihilates electric 
field quanta. There is, however, no utility in explicitly 
introducing E and A operators, at least at the mean field 
level we proceed with in the following, and we will not 
do so. 

It is convenient to introduce a rotor variable ip r canon- 
ical conjugate to Q r ; 

[<Pr,Qr]=i, (9) 

which gives 

$ r = e - <Vr , $t$ r = l. (10) 

Using the representation in Eqs. ([5])- Q, we rewrite spin 
Hamiltonian Eq. 



r r p,^jtv 
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+ i? r e„ s r!r+?7 r e l , + fl.C.J +COnst.. (11) 



The total "electric charge" Q = J2 r Qr commutes with under the gauge transformation, 
the Hamiltonian Hqed- Furthermore, Hqed is invariant 



(12) 
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This completes the reformulation as a compact U(l) 
gauge theory with bosonic matter: a compact interacting 
Abelian Higgs model. 



III. MEAN-FIELD THEORY 



distinguish this from ordinary Curie- Weiss mean field 
theory, we denote this treatment as gauge Mean Field 
Theory (gMFT).ES Specifically, we decouple the various 
terms in Eq. |TT| ) as follows: 



We now proceed with a mean- field analysis of the gauge 
theory in Eq. (11), following and extending Ref. IT51 To 
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The second decoupling (Eq. ( |T4| )) is introduced here (and 
was not considered in Ref. I15p to deal with the J±± inter- 
action, which involves not only interaction between the 
"gauge fields" (s r)1 ./) and spinons, but also between the 
spinons themselves (being quartic in $ r , $J operators). 

After the above decoupling, some ansatz must be made 
to determine the form of various expectation values. Of 
particular interest are the "magnetic" expectation values 



= l<s 



± ,) 

r.r' I 



(15) 



which are formally similar to the bond operators appear- 
ing in large N slave particle theories of quantum antifer- 
romagnets. In particular, the phase A r yof this expecta- 
tion value has the physical interpretation of a sort of "av- 
erage" gauge field experienced by the spinons. Though 
this is not itself gauge invariant (see Eq. (12)), the net 



flux of this gauge field, or equivalently the product of 
(s+ r ,) expectation values around any closed loop, is phys- 
ically meaningful. Different patterns of this flux describe 
different QSL states, with different Projected Symmetry 
Groups, or PSGs.^S Within the gMFT formalism, the 
choice between different flux patterns, or PSGs, should 
be made based on comparison of the mean-field free en- 
ergy for different Ans&tze. 

Instead, we choose the flux pattern based on a non- 
mean-field, but perturbative argument. This has the ad- 
vantage of being simpler, and also correct beyond mean 
field theory in the perturbative regime, but could in prin- 
ciple break down by missing some other phase at larger 
coupling. We leave a more exhaustive study of energetics 
of different PSGs as an open problem for the future, but 
expect that the choice made here is probably correct in 
most cases of interest. 



In the perturbative regime, J±±/J zz < J±/J zz <c 
1, the leading contribution in degenerate perturbation 
theory^ occurs at order (J±/J zz ) (the contribution from 
J ±± is subdominant), and gives an term in the Hamil- 
tonian proportional to the sum of the cosine of the flux 
through each hexagonal plaquette of the dual diamond 
lattice, 



J^cos(V x A), 

O 



(16) 



where (V x A)q = ^2 ie Q A TuTi+1 denotes a lattice 
curl of the gauge field A TT >, i.e. the flux through the 
hexagon containing the sites i. This calculation deter- 
mines the PSG to be: (1) for J± > 0, a 0-flux state, with 
cos(V x A) = 1 or (2) for J± < 0, a 7r-flux state with 
cos(V x A) = —1. This is somewhat consistent with the 
physical intuition that for J± < 0, the XY pseudospin 
order is frustrated, requiring formation of a more com- 
plex ground state. Indeed, we will see later that quantum 
fluctuations are greatly enhanced in 7r-flux state, leading 
as a consequence to much enhanced domain of stability 
of the QSL phase relative to the case J± > 0. In the 
following, we will consider the zero-flux and 7r-flux cases 
separately. Fig(2] shows mean-field ansatz of gauge fields 
for zero-flux and 7r-flux cases. Thick links have A r r i = tt 
and other links have A r r > = 0. There exist many other 
gauge field configurations that produce the same physical 
state. 
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FIG. 2: Mean-field ansatz for gauge fields that preserve the 
symmetry of Hamiltonian Eq. (Ill: (a) 0-flux state, and (b) 



7r-flux state. In (b), the thick links have A r y — n and other 
links have A r y = 0. Thick solid (dashed) links are aligned in 
upper (lower) hexagonal plane perpendicular to (111) direc- 
tion in a diamond lattice. 



IV. FM CASE (J± > 0) : ZERO-FLUX STATE 

In this section, we carry out an analysis of the "fer- 
romagnetic" case, J± > 0, for which the 0-flux state is 
stabilized. The zero flux state was considered in Ref. IT51 
but with J±± = 0. 



A. Order parameters 

The gMFT treatment introduces several self- 
consistently determined "order parameters" , whose 
interpretation requires additional care in comparison to 
ordinary mean field theory, owing to gauge redundancy. 
We discuss this now. In the "gauge" sector, there are 
two types of expectation values: (sj p ,), which is gauge 
invariant and direct ly p roportional to the local magnetic 
moment, see Eqs. (2 6), and (s„;), which as discussed 
above is not gauge invariant and related to the "average" 
and fluctuations of the magnetic vector potential. There 
are several other order parameters in the "matter" 
sector. While it is not explicit as a decoupling in the 
gMFT scheme, the spinon condensate itself, ($ r ), is 
an important (non-gauge invariant) order parameter. 
In comparison to the prior case, the decomposition in 



Eq. (14) introduces two new types of order parameters. 
First, there are pairing terms between two spinon 
fields on the same sublattice, ($ r+rtre ^& r+rirev ). This 
composite held carries electric gauge charge 2. Second, 
there are particle-hole terms which are gauge neutral 
but connect the two sublattices, ($J$ r +^ re ). The 
possibility of non-zero expectation values of these fields 
enlarges the spectrum of phases which may occur within 
gMFT. 

For the FM case, we have found that the mean field so- 
lutions do not break translational symmetry. Presuming 
this to be true, we can classify the different possible types 
of solutions by their order parameter expectation values, 



and we discuss the physical meaning of these patterns 
now. The different possibilities are listed in Table. |LT| 
and summarized below. 



1. Coulomb phases 

Within a mean field description, Coulomb phases are 
those in which the U(l) gauge symmetry is unbroken by 
charged condensates, i.e. ($ r ) — ($ r+)7i . eti $ r+ , ;r e„) = 0, 
and where the gauge fluctuations are not so strong as 
to wash out the average of the gauge elements, i.e. 

( s r,bfr') # °- This leaves ( S ?,r') and (^r^r+^e,,) un- 
determined, and depending upon their values, different 
phases can be realized. 

a. U(l) QSL - When (s* r ,) = ($ r $ r+ „ re J = 0, 
all physical (global) symmetries of the system are unbro- 
ken. Due to the non-zero value of (s rr ,), the spinons 
are able to coherently hop and propagate. Thus this 
phase may be characterized as a U(l) QSL with propa- 
gating spinons and an emergent U(l) gauge field (which 
appears in the mean field treatment by fluctuations of 
the phase A r y). This is the realization in gMFT of the 
QSL discussed perturbatively in Ref. [3 

b. Coulombic ordered phase When either 
($ r $ r ') 7^ 0, (s* r ,) 7^ 0, or both, global physical 
symmetries of the pscudospin Hamiltonian are broken. 
For instance, (s P r , ) ^ implies time- reversal symmetry 
breaking and the presence of spontaneous magnetic 
dipole moments oriented along the local (111) axes. 
However, phases of this nature are not trivial ordered 
states, since like the QSL, they host propagating decon- 
fined spinons and an emergent gapless Coulomb gauge 
field. A particular example of this type, dubbed the 
"Coulomb ferromagnet" , was discussed in Ref. HH but 
more general such phases may occur when the full phase 
space is considered. They do not, however, occur in the 
gMFT solution in the following subsection. 



2. Confined phase — Ising order 

Another possible phase is one in which the gauge fluc- 
tuations are sufficiently strong that they destroy the bond 
operator expectation values, (s*,.,) = 0. In this case, the 
amplitude for spinons to hop vanishes, and so they do not 
coherently propagate. This should be considered there- 
fore a confined phase. Since the s r r ' are spin-1/2 oper- 
ators, the only state consistent with the above condition 
is polarized along the z direction, so that (s pr ,) ^ 0, im- 
plying broken time-reversal symmetry and Ising magnetic 
order. In this state, the gauge fluctuations are insignifi- 
cant, and there are no emergent low-energy gauge fields. 
Indeed, such a state would be considered a classically 
ordered spin ic d 33 * 34 [ In our calculations, we find that 
this phase does not occur in the region of phase space we 
studied. 
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3. Higgs phase — quadrupolar order 

It is well known that in gauge theories there are two 
ways to remove the gauge fields from the low energy 
physics: confinement, as described above, and the Higgs 
mechanism. In the Higgs mechanism, a boson carrying 
the fundamental gauge charge condenses, leading to a 
"Meissner effect" for the gauge flux and a gap for the 
photon. Though the Higgs mechanism for the transition 
seems very different from that occuring in confinement, 
it is believed that, in the absence of global symmetries, 
the confinement and Higgs phases are indistinct, and can 
be adiabatically transformed into one another. There- 
fore, like the confined state discussed above, the Higgs 
phase (s) which occur in gMFT are conventional states 
of matter, absence any exotic excitations or topological 
properties. In our case, however, the Higgs phases can 
be distinguished by symmetry from the confined one. 

To see this, consider a state with spinon condensation 
at wavevector k, i.e., 

(*r) = (<J>k)e lk ' r ^ 0. (17) 

For the spinons to condense, they must obviously propa- 
gate, so we must also assume (s^ r ,) ^ 0. The mean field 
result is then long-range ordering of pseudospins on the 
XY plane (which corresponds to quadrupolar ordering in 
physical terms) since 

(S + ) = (d>t)(s+ r+ e (1 >( ( I ) r + e (1 >^0. (18) 

We indeed find Higgs phases of this type in the solution 
of the gMFT equations given below. 



4- Charge 2 Higgs phase — Z2 QSL 

The remaining possible phase is one in which the 
charge 2 Higgs (spinon pair) condensate is non-vanishing, 
($ r +, reji <& r +ri r e v ) 7^ 0, but the fundamental spinon is un- 
condensed, ($ r ) = 0. In this case, the U(l) gauge sym- 
metry is broken down to Z%, and the resulting state is a 
gapped QSL with only Ising-like gauge charges. Spinons 
remain deconfined, but become mixtures of particles and 
anti-particles, much like Bogoliubov quasiparticles in a 
superconductor are mixtures of electrons and holes. Due 
to the absence of spinon condensation, this state generi- 
cally has vanishing spin expectation values and need not 
break symmetries. Like the confined and Coulombic or- 
dered phases, this state does not, however, occur as a 
ground state in gMFT, in the parameter regime we have 
so far studied. 



B. Self-consistent equations 
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TABLE II: Classification of possible phases occuring in the 
gMFT treatment. 



where H s describes the gauge sector as decoupled "spins" 
s r r ', and the spinon part £/$ contains only <fr r ,$J and 
Q r operators and is quadratic in them. The gauge part 
H s is trivially soluble, as different bonds are decoupled, 
and each s rr < is then treated simply as a s = 1/2 spin in a 
field. The spinon part, however, is non-trivial, and actu- 
ally still strongly interacting, since $ r is actually defined 
in terms of the fundamental rotor field ip r , c.f. Eq. (10 1. 



Once the replacements in Eqs. ( 13|14 ) have been made, 
the mean field Hamiltonian reduces to a sum of Hamil- 
tonians, 

H —> H g MFT = H s + (19) 



To proceed with it, we follow Refs. 36 15, and replace 
the "hard" constraint on $J$ r = 1 with a "soft" an av- 
erage one, enforced by a Lagrange multiplier A r . This 
is equivalent, quantum mechanically, to promoting the 
spinon field to a complex rotor, with <I> r = x r + iy r and 
Qr = Px r + Wy r , where x and p variables are canonically 
conjugate coordinates and momenta. After this substitu- 
tion, ff$ becomes quadratic and soluble. The Lagrange 
multiplier, appearing as a mass for $ r , is adjusted to 
maintain (<f>J$ r ) = 1 on every site. 

In the following, we use this formulation to calcu- 
late the necessary expectation values and impose self- 
consistency. We assume here a zero flux state, and no 
breaking of translational symmetry. We also neglect the 
J z ± coupling, in which case one can show that the en- 
ergy is minimized when (s^ r ,) =0. Then in general there 
are many variables: 4 distinct gauge fields (s+ r ,) on the 
four orientations of diamond bonds, 8 spinon pair fields 
(2 on a single site and 6 distinct orientations of pairs 
connecting same-sublattice sites in different unit cells), 
4 A-B sublattice mixing field on diamond bonds, and 2 
Lagrange multipliers, one for each of the two basis sites. 
This makes 18 distinct mean field parameters, a general 
analysis of which is daunting. To proceed, we looked for 
self-consistent solutions with fewer parameters. We em- 
ployed a rather general ansatz, containing both pairing 
and A-B sublattice mixing, but imposing some discrete 
symmetry constraints. We discuss the comparison of the 
energy of these solutions in the subsequent subsection. 
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1. Mean-Field Ansatz 

We introduce the following ansatz including both pair- 
ing and A-B sublattice mixing: 



A = 


( S r,r±e M )' 




(20) 








(21) 


Xo S = 


($r B $r B >, 
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xt = 


($r B -e tl $r B -eJ, 


e M — e„ € j'fe plane, 


(23) 


xf = 




e M — e„ € j'fc plane, 


(24) 


Cm = 






(25) 



As mentioned in the previous section, the gauge sector is 
trivially soluble, leading to A = 1/2. The spinon action 
part can be rewritten in matrix notation: 



duj r< 



2tt ^ K V 

k>0 v 



2J Z 



(26) 



where 



( ^ \ 



4>u 



M(k) 



<J>k 



/ A u (k) A 12 (k) 

AJa(k) A n (k) 
C*(k) 

V 



C(k) 

C*(-k) 
B u (k) S 12 (k) 



(27) 



(28) 



C(-k) B? 9 (k) Bn(-k)7 
and A a ^(k), B a ^(k) and C(k) are defined as 
An(k) = B n (k)= - r \- \ ' P -* k < e 



- J±A 2 £ e 



A i2 (k) 
C(k) 



J±±A 2 
2 

J±±A 2 
2 

J±±A 2 



ni/X 



/ 1 ' ' 



- 7n„Xo e 



-ik-(e t , 



M#^ 



ik-e„ 



(29) 



Finally, to render the mean field problem solvable, we 
replace the constraint |$ r | = 1 by the "softened" con- 
straint -k ^2 r \& r \ 2 = 1, and implemented the latter by 
including a Lagrange multiplier term into the action 5$ . 

Using this formulation, the mean field Hamiltonian al- 
lows one to calculate (Hqed) (Eq.(ll)) and minimize 
this variational energy. We found and compared several 
self-consistent solutions of the gMFT equations, which 
are subsets of the general ansatz given above. First, we 
considered two limits allowing for pairing, or A-B sublat- 
tice mixing, but not both: 



(i) 
(ii) 



6. = o, 



xT ] ¥> o, 



yA{B) 



MB) 

X 



Y A(B) 



(30) 
(31) 



While self-consistent solutions may be found for both 
these cases, we find that the minimum energy solutions 
always have either vanishing pairing/sublattice mixing 
(i.e. describe the U{\) QSL) or exhibit spinon condensa- 
tion. 

However, the both condensed solutions are unnatural, 
insofar as once a single $ field is condensed, all the ex- 
pectation values Xo^ B i xt^ B > Cm would be expected to be 
non-zero. Guided by the above cases, we found a self- 
consistent ansatz where all these were allowed to be non- 
vanishing, with the relations \o = Xo j E^T/^X^L = 
J2^ v %vX% ± and f = 6 = ~ii = ~ih 0, for 
{i,j, k} and permutation of {1,2,3}. This more general 
ansatz describes both condensed and uncondensed states, 
and was found to capture all the physical minimum en- 
ergy solutions. 



2. Spinon condensation 

In the gMFT scheme used here, Higgs phases in which 
the single spinon field is condensed, (<& r ) ^ 0, also occur. 
This may appear surprising since the single spinon field 
was not introduced explicitly as an order parameter - see 
Eqs. (13) and (14). Instead, spinon condensation occurs, 
as discussed in Ref. 15, via the same mechanism as does 



Bose-Einstein condensation in the non-interacting Bose 
gas. In particular, when a condensate is present, the La- 
grange multiplier A adjusts itself self-consistently so that 
the minimum energy spinon state lies, in the thermody- 
namic limit, at precisely zero energy. For large but finite 
volume, a non-intensive part of the A leads to and con- 
trols the condensate, manifesting itself via off-diagonal 
long range order in the spinon Green's function. This is 
discussed in more detail in Appendix | A 1 b Captured in 
this way, spinon condensation does not introduce any ad- 
ditional self-consistent variables, and only requires care- 
ful treatment of any zero energy modes and the infinite 
volume limit. This in turn means that the above ansatze 
describe Higgs phases as well, for appropriate values of 
parameters. 



C. gMFT Phase Diagram 

We minimized the variational energy using the above 

for the formula- 



re self-consistent 



ansatz numerically (see Appendix Alb 
tion of the variational energy) . In fact 
gMFT equations are solved for any local minima of the 
variational energy, so it is sufficient to search for the 
global minimum of the latter. That determines the T = 
phase diagram as a function of J±/ J zz > and J±±/J zz 
(we assume J zz > throughout). Note that by a canon- 
ical transformation, S 1 * 1 — > ztiS ± , we can always choose 
J±± > 0, without loss of generality. The results are 
shown in FigjSJ 

The full phase diagram contains three distinct phases 





0.8 




0.6 






+i 
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0.2 



noncoplanar FQ 


U(1)QSL 
Spin Ice 


AFQ 
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FIG. 3: Phase diagram of two dimensionless parameters 
J±/Jzz vs J±±/ J zz . Four distinct phases exist : classical spin 
ice (at the origin), U(l) QSL, AFQ and FQ. (more details in 
the main context) 



in addition to the classical point corresponding to 
the nearest-neighbor spin ice: a deconfined U(l) QSL 
phase and two Higgs phases, corresponding to XY 
ferro-pseudospin (antiferro-quadrupolar) and antiferro- 
pseudospin (noncoplanar ferro-quadrupolar) orders. Un- 
fortunately, the Z2 spin liquid phase with non-zero pair- 
ing but a spinon gap is never the minimum energy solu- 
tion. The QSL or Coulomb phase occurs in the small 
J±,J±± region, consistent with perturbative expecta- 
tions. In this model, infinitesimal J± and/or J±± interac- 
tions "melt" the classical spin ice, creating a dynamical 
"photon" excitation and emergent quantum electrody- 
namics. This phase is found to be more stable against 
J_l__l_ than to J z ± , the latter having been studied already 
in Ref . QH 

The Higgs or ordered phases merit some further de- 
scription. With increasing J±/J zz but J±± — 0, the 
U(l) QSL phase remains stable untill j^| c ~ 0.19, at 
which spinons start to condense at a wave vector ko = 
for both A and B sublattices. This induces a classical XY 
order categorized in Table|n]and has the ordering struc- 
ture shown in Figj4] (a). This phase has already been 
obtained by a classical MF analysis^, and in gMFT for 
J±± = (P- From Eqs. (fTf} and (pi, the spinon conden- 



sate at ko yields a ferroic ordering of the XY component 
of pseudospins, for instance, given by 



(32) 



for pseudo-spin on sublattice i. It spontaneously 
breaks the threefold rotational symmetry while the 
twofold rotational symmetries are preserved. This 
ferro-pseudospin ordering structure is interpreted as an 
antiferro-quadrupolar order for Pr 3+ case as is clear from 

i = 0. Namely, it 




FIG. 4: (a) AFQ : an /-electron charge distribution under 
the antiferro-quadrupolar order induced by a spinon conden- 
sation ($) 7^ at k = (000) wave vector. Here, we have 
taken the Pr 3+ case with the local ground-state non-Kramers 
doublet a\4a) - f3a\a) + 7] - 2a) with a « 0.970, /3 « 0.075 
and 7 « 0.230, where \ J Z ) is an eigenstate of the z compo- 
nent of the total angular momentum in the local framfP. (b) 
Noncoplanar FQ : an /-electron charge distribution under the 
noncoplanar ferroquadrupolar order induced by a spinon con- 
densation ($) / at k = 2tt(100) for the Pr 3+ case, (c) 
Antiferromagnet : an XY antiferromagnetic ordering induced 
by a spinon condensation {<&) 7^ at k = (000) wave vector 
in the case of half-integer spins, (d) Noncoplanar ferromag- 
net : an XY noncoplanar ferromagnetic ordering induced by 
a spinon condensation ($) 7^ at k = 27r(100) in the case of 
half-integer spins. 



produces an /-electron distribution shown in Fig[4](a). 
When J±± > is sufficiently large and J± is small, the 
QSL becomes unstable to a different Higgs phase, with 
spinon condensation at ko = 27r(100) or the symmetry 
related points, on both A and B sublattices. Note that 
quantitatively the QSL phase is wider in the J±± direc- 
tion than in the J± one: ^f^\c ~ 0.31, compared to 

7^|c ~ 0.19. This suggests that the U(l) QSL phase is 
more stable against the J±± interaction than the J± in- 
teraction. This can be understood from degenerate per- 
turbation theory. The J±± interaction induces a non- 
trivial contribution only at the sixth order, 0(J±±/ J z ), 
whereas the comparable term is induced already at the 
third order in J±. As for the other Higgs phase, the or- 
dering structure is understood again from Eqs. (17) and 
( 18 ). One of the symmetry-broken ground states is 



Eq. ([3| and the relation ^ 



((So), (§!>, (S 2 ), (S 3 )) « l^fifayu-m, -y 3 ). (33) 

It spontaneously breaks both the threefold rotational 
symmetry and the cubic symmetry, and loses two of the 
twofold rotational axes. This antiferro-pseudospin struc- 
ture is interpreted as a noncoplanar ferro-quadrupolar 
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order in the Pr 3+ situation, as is clear from Eq. ^ and 
the relation (y + yi — y 2 — §3) | (100). It creates an 
/-electron distribution shown in Figj4] (b). It is worth 
to note that the above ferrro- and antiferro-pseudospin 
structures, associated with the antiferro- and noncopla- 
nar ferro-quadrupole orders shown in Fig[4](a) and (b) for 
Pr 3+ cases, are directly related to XY-magnetic orderings 
shown in Fig. 4 (c) and (d) when we consider half-integer 
spin rare-earth pyrochlores instead of integer spin case. 



D. Phase transitions 

Within gMFT, the phase transition between the U(l) 
QSL and AFQ state is second order, as indicated by 
a continuous change of the MF variables \v from zero 
to finite values across the phase boundary (solid line in 
Figj3}. A low energy continuum action for this transition 
is simply an Abelian Higgs theory, with a charged bosonic 
matter field (representing the condensing spinons) cou- 
pled to a dynamical gauge field A. When gauge fluctua- 
tions beyond the mean field are included, such transitions 
are usually driven weakly first order: 35 It is interesting 
that, within gMFT, the phase boundary between QSL 
and AFQ phases is precisely vertical, as seen in Figj3] 
This is because in the QSL phase arbitrarily close to the 
phase boundary, both spinon pairing and A-B sublattice 
mixing is absent, so that the J±± interaction gives zero 
contribution to the energy. 

By contrast with the above case, we find that the QSL 
to FQ transition is strongly first order already in gMFT. 



This is indicated by the dotted line in Fig(3j Fluctua- 
tion effects will not change this conclusion. Note that 
this phase boundary has a positive slope, i.e. the FQ 
state is suppressed by increasing J± . This is because J± 
interaction prefers instead the AFQ state. 



V. AF CASE (J± < 0) : tt-FLUX STATE 

A similar analysis can in principle be completed for 
the case J± < 0, for which the XY-pseudospin order 
is frustrated. This case, however, introduces significant 
new complexities which are beyond the scope of the 
present work, and will be discussed in a future publi- 
cation. Here, we confine ourself to the line in the phase 
diagram J±± = in the AF pseudospin region. 

As discussed in Sec |IV| the J± < favors a 7r-flux 
state, in which all the hexagons carry a flux V x A = tt 
(mod 27t). For calculations, it is necessary to choose a 
gauge with a specific assignment of A r r i having it flux, as 
shown in Fig(2| (b). In this case, the unit cell is doubled 
compared to the case of FM J± and contains four sublat- 
tices (comprising 2 diamond sites in each of the 2 mag- 
netic unit cells ). This gauge field pattern can be repre- 
sented as A Y , x+e ^ = e M Q ■ r where (e , ei, e 2 , e 3 ) = (0110) 
and Q = 27r(100). In the QSL state, this leads to 
( s ^ r +e fl ) = Ae"^- r , with A = 1/2. 

In this fixed gauge, we consider the spinon dispersions 
for J±± = 0. Within gMFT (see Eq. Q), the A and B 
sublattices are decoupled and the spinon action is 



S$ = 



f dT E ^ T $;a$ r + A A ^(|$ r | 2 -i) + A B ^(|$ r | 2 -i) 

^ r<£A,B z r<£A reB 

+ J ±{ J2 J2^r+^+^Kr+ e J(4,r+ e J + E ( # S-e„ ^-e, ) (*r,r-e„ ) } 

r£A ft^v r£B fx^v 

d 3 k dw 



V BZ 2tt 



( aft 



aft 



2J« 



X A + P A -P. 2 A - iPA 



Pf + iP^ 7TT \~ \ A — Pi J \ $ A 



(A^B) 



(34) 
(35) 



3 2J„ 
I 



k,2 



Here 



Pi = 4J±A 2 cos cos 



P A = 4J±A 2 sin 



k,. 



sm 



2 2 

?A at A 2 . kr. . K 



P3 = 4J±A^ cos — sin 



2 ' 



pB 



(36) 



with P/ 3 = Pf(k -> k + 7r(lll)), 
k + 7r(lll)) and P 3 B = P 3 A (k -> k + Tr(lll)). We now 
seek a QCP between the U(l) 7r-flux QSL and a magnet- 
ically ordered phase. Since the A and B sublattices are 



decoupled, it is sufficient to focus on one sublattice, for 
instance, A. As before, we adopt a "softened" constraint 
lEri^'^r 1 ) = 1j which leads, assuming no spinon 
condensation, to 



IJxz 1 f d 3 k 
2J±AJ V BZ 



= 1. 



2A + 2JA 2 



(37) 



PI 



Pi - 



Here we defined A = \ A / J±A 2 and 
S a =i 2 3 \_Pak/J± A 2 ] 2 ■ The spinon condensation point 
occurs when the integrand diverges, which gives A c = 



10 



-y/maxk-P^ = 4. By substituting A c = 4 and evaluating 



the integration in Eq. (35) at this point, we obtain 



\J±\ 

J 7.7. 



4.13. 



(38) 



This is the main result of this Section. We observe that 
the QSL phase is much more stable to antiferromagnet 
J± than to ferromagnetic J± . This is rather natural since 
the competing XY pseudospin order is frustrated in the 
antifcrromagnetic case. We can understand this more an- 
alytically from the spinon dispersion in the 7r-flux state, 
which has the form E£ = ^(P^) 2 + (P^) 2 + (^ 3 k) 2 - 
This form, which describes states in either A and B sub- 
lattice, has a degenerate set of energy minima, consist- 
ing of lines in reciprocal space (e.g. E£ is minimized 
for k = (A;, 0,0) with an arbitrary real number k, and 
there are several other similar minimum energy lines). 
In contrast, in the FM case, k = (000) uniquely gives 
the minimum energy. This effectively lowers the spatial 
dimensionality at low energies, increasing the stability of 
the U(l) QSL. We note, however, that this line degener- 
acy is emergent and is not protected by any symmetry. 
Effects beyond gMFT should be taken into account to 
further split this degeneracy. Such effects would be es- 
sential in determining the nature of quadrupolar ordering 
in the Higgs phase beyond the critical point. We expect 
this physics to lead to a significantly richer phase diagram 
when J±± interaction is included. 



VI. SUMMARY 

In this paper, we have studied the generic pseudospin- 
1/2 model describing nearest-neighbor coupling of 
ground state magnetic doublets of rare earth ions on the 
pyrochlore lattice. We showed how to extend the gMFT 
treatment of Ref. [TS] to take into account all the sym- 
metry allowed interactions, which requires a significant 
extension for the method. We focused on the case of a 
non-Kramers ion, for which three interactions exist: an 
Ising spin-ice interaction J zz (which we presume always 
takes the non-trivial frustrated sign), a U(l) symmetric 
planar exchange J± and an in-plane anisotropic exchange 
J±± ■ For the case of ferromagnetic symmetric planar ex- 
change, we obtained a complete gMFT solution. This 
situation favors "zero flux" states in the gauge theory 
formulation. We obtained a finite region in the phase 
diagram supporting a U(l) quantum spin liquid (QSL) 
state, described as a type of emergent quantum electro- 



dynamics. Phase transitions from the U(l) QSL state to 
two types of planar pseudospin orders were found to oc- 
cur by the Higgs mechanism with increasing J± and J±± ■ 
For large J± this yields an antiferro-quadrupolar phase, 
while large J±± yields a ferro-quadrupolar state. In the 
case of antiferromagnetic symmetric planar exchange, a 
7r-flux state is preferred in the gauge theory, and the gen- 
eral solution was too complex to attempt here. However, 
we did prove that the increased frustration in this regime 
greatly increases the stability of the U(l) QSL state. 

It is hoped that these results form some basis for under- 
standing experiments in the non-Kramer's pyrochlores 
Pr 2 TAf 2 7 with TM =Sn, Ir, and Zr. Future studies 
should complete the full phase diagram in the antiferro- 
magnetic planar symmetric exchange case, and address 
the accuracy of the gMFT results by comparison with 
other methods, considering the role of further neighbor 
interactions, lattice distortions, and disorder. 
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Appendix A: Mean-field approximation 

In this section, we proceed two steps of mean-field 
(MF) approximation to make Eq. (Ill soluble. 



1. Green's functions and energy 

a. Green's function 



Using the MF ansatz listed in Eqs. (20)-(25), the 



Hamiltonian of particle part are written as 
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H„ — 



H„ — 



J 



2 Y.Ql + H P 

r 



J±±A 2 



r I 



E { E ^ ) $r +f)r e t $r + , )r e„ + E IJ^X^&M + h.C 



-2J±± A 2 ]T ^ X! 7;>^^^ r+J?r e„ + ft.c 

r I 



(Al) 



(A2) 



with r + = B and T_ = A, which leads to the action, 
S p = I dr ]T -^9 T ^G> r $ r + H. p 

J r£A,B zz 

+ A r (l*r| 2 -1)) (A3) 

reA.B 



In Eq. ( A3 1 , the first term comes from integrating out Q T 



and the last term is for Lagrange multiplier which con- 
strains |$ r | 2 = 1. In large N limit, the integrals become 
sharply peaked at the saddle point, say X A< - B \ Hence 
we pull out A,, from the summation with its saddle point 
value \ A ( B \ This is consistent with softening local con- 
straint |$| 2 = 1 to its average J2 r l^ r| 2 = N. Using this 
saddle point approximation, Eq. ( |A3[ ) can be rewritten in 
a Fourier transform and this results in Eq. ( 26 ) . Fourier 
transform of (f> r T is defined, 



1 



Nu.cJ 2w 



E* 



-i{u)T— k-r) 



(A4) 



where k is wave vector, uj is an imaginary frequency and 
N u _ c is the number of unit cell. Then Green's functions 



are represented as 



uj 2 + 2J zz (X + e m ) 



(A5) 



The Matsubara sum of frequency leads 
f duj 

g af 3(k) = / — g a p(k,u) 



E 



2J zz <j>™*<P™ 



2y / 2J zz (X + e m ) ^ u„ 



(A6) 



where <j)™ is the a th component of m th eigenvec- 
tors for M, e m is t he m th eigenvalues for M and 
^2J ZZ {\ + 

b. Variational energy 

We consider the energy {Hqed) by taking an expec- 
tation value of Eq. ( 11 ). 



(Hqed) = J -f E<Q 2 > - J± EE(C,^><s,"X.J<sXJ 

r r fi^jtii 

(A7) 

I 



First of all, let's consider J zz /2 ^2 r (Qr)- As we men- 
tioned in Sec jlV B| this term can be represented as 
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/ ^E( 2 -^ 2 (^( k ' w ) + ^3(k,c)) 

k 



k m 



Here, we used (p 2 ) =1/ZJ dpdxp 2 e - 1 dr(p 2 +i P ±+f( x )) 
1 — (i 2 ) where £T is partition function. 



Finally, Eq. (A7| can be rewritten as, 



(A8) 
(A9) 
(A10) 
(All) 



(Hqed) = EEj' 



k rn 



' 2JV„. r 



J±±A 



k^ tri' 



2 I 

) + E E ^-f + A U ro 

km m 

[{EE E E E ^e—-'-.- 



{ E E U^r^ E t,-**^- } { E E Itt^'c' } 



2J Z 



kA m 



k^ 77l' 

2J Z 



2A ^ [ E *v { E E £ «■ ) { E E £ c'c'e"— } 

li^v V.a m k^ m' 



C.C 



(A12) 



2. spinon condensation 

When spinons condense at momentum ko, summation 
of k in the first Brillouin zone can be replaced by 



-E2« 



g(fco) 



- E f( k ) ( A13 ) 



k#k„ 



Spinon condensation also affects to a lagrangian multi- 
plier term and leads A to be 



A = A + 



A' 



N, 



(A14) 



where Aq is the minimum of e„ 



N, 



^^/(k,, m (A,k))c(k)^(k)4 /(k ,^(A',k ))c*(ko)^ l (ko) 

d 3 k 



k m 



/ ^£/(k, Wm (A„,k))c*W(k) (A15) 



m is the m th eigenvectors of M which has the minimum 
eigenvalue of e m . 
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